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Abstract 

It was recently observed that sand flowing down a vertical tube sometimes 

forms a traveling density pattern in which a number of regions with high 

density are separated from each other by regions of low density. In this 

work, we consider this behavior from the point of view of kinetic wave theory. 

Similar density patterns are found in molecular dynamic simulations of the 

system, and a well defined relationship is observed between local flux and local 

density - a strong indicator of the presence of kinetic waves. The equations 

of motion for this system are also presented, and they allow kinetic wave 

solutions. Finally, the pattern formation process is investigated using a simple 

model of interacting kinetic waves. 
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Systems of granular particles (e.g. sand) exhibit many interesting phenomena, such as 
segregation under vibration or shear, density waves in the outflow through tubes and hop- 
pers, and probably most strikingly, the formation of heap and convection cell under vibration 
In granular flows through a narrow vertical tube, Poschel found that the particles 
do not flow uniformly, but form high density regions which travel as coherent structures 
with a velocity different from the center of mass velocity. He also reproduced these density 
waves using molecular dynamics (MD) simulations f^. However, the motion of these high 
density regions and the mechanism which is responsible for their formation are not fully 
understood. 

In this paper, we present numerical and theoretical evidence that these density waves are 
of a kinetic nature 0. Using MD simulations, we measure the dependence of the particle 
flux on the density. We find a well-defined flux-density relation - an indication that a 
kinetic wave theory describes the behavior. A direct measurement of the velocity of these 
high density regions shows a dependence on the mean density which is in good agreement 
with the predictions from kinetic wave theory. On the theoretical side, we consider one 
dimensional equations of motion for the density and the velocity fields in the tube. These 
equations, together with Bagnold's law for friction |]^, allow kinetic density wave solutions. 

In order to understand the formation of these high density regions, we consider the 
general problem of interacting kinetic waves. We first show numerically that a system with 
an initially random density field evolves to a configuration in which neighboring regions have 
a high density contrast. At the early stage of development, we can show analytically that 
the density contrast between nearby regions increases linearly with time. 

We first discuss the MD simulations of the system, and begin with a brief description of 
the interparticle force laws that were used in our calculations. The particles interact with 
each other (or with a wall) only if they are in contact. The force that acts on particle i 
due to particle j can be divided into two components. The first, -Fjl^,, is parallel to the 
vector r = Ri — Rj, where Ri and Rj are the coordinates of the centers of particles i and j 
respectively. We refer to this as the normal component. The second component, orthogonal 
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to r, is the shear component F^^^. The normal component is given by 

F^^i = K{ai + aj - |r|)3/2 _ ^^^^^LL^ (la) 

where ai{aj) is the radius, and mi {rrij) the mass of particle i (j). Also, rrie is the effective 
mass rriimj/ {mi + mj), and v = df/dt. The first term in Eq. ( [Ta| ) is the Hertzian elastic force, 
where kn is a material dependent elastic constant. The second term is a velocity dependent 
friction term, where 7^ is a normal damping coefficient. The shear component is given as 

= -7sme^, (lb) 

where s is defined by rotating r clockwise by 7r/2. The shear force, Eq. (ITbl), is simply 
a velocity dependent friction term similar to the second term in the normal component. 
Finally, we must specify the interaction between a particle and a wall. The force on particle 
i, in contact with a wall, is given by Eqs. (1) with aj = 00 and rrij = 00. The choice of the 
interactions defined by Eqs. (1) is rather typical in the MD simulations of granular material 
P]. A detailed explanation of the interaction is given elsewhere Q. 

For simplicity, we study granular flows in 2 dimensions and use a fifth order predictor- 
corrector scheme to integrate the equations of motion, calculating both the positions and 
velocity of each particle at all times. The tube is modeled by two vertical sidewalls of length 
L with a separation W, and we apply a periodic boundary condition in the vertical direction. 
Between the sidewalls, particles of radii 0.1 are initially filled with a uniform density of po 
(throughout this paper, numerical values are given in CGS unit). The particles begin to 
move under the influence of gravity, and soon reach a steady state, where the gravitational 
force is balanced by the frictional force from the interactions with the sidewalls. 

In Fig. |1|, we show the time evolution of the density and the velocity fields for L = 15 
and W = 1, measured at every 5 ms. At a given time, we divide the tube into 15 vertical 
regions of equal length, and measure the density and the average velocity in each region. 
These fields are displayed as a vertical strip of square boxes, where each box corresponds to 
a region in the tube. The grayscale of the box is proportional to the value of the field in that 
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region. The parameters we used in this simulation are kn = 1.0 x 10^, 7„ = 7s = 5.0 x 10^, 
with the time step 5.0 x 10^^. The initial density po is 25 particles per unit area. In the 
figure, we find (1) a region of large density fiuctuations is formed out of the initially uniform 
system, (2) the fiuctuations seem to travel with almost constant velocity (different from the 
center of mass velocity), and (3) there seems to be strong correlation between the density 
and the velocity fields. These findings remain true for the simulations we have performed 
with different values of 7, kn and po, except when po is very small, where a steady state 
is not reached. These traveling density patterns were first observed in the simulations by 
Poschel §. 

In order to quantitatively study the correlation between the density and the velocity 
fields, we measure the local particle fiux as a function of the local density in the following 
manner. Once the system has reached a steady state, we measure the mean velocity Vi and 
the density pi in region i. The fiux j{p) for a given density p is then taken to be p ■ {v{p)), 
where () is a time average over all regions which had a particular density p. The fiux- 
density curve, obtained by averaging over 10,000 iterations, are shown in Fig. 0. Here, the 
parameters are the same as those of Fig. |I|. The fact that a well-defined fiux-density curve 
exists suggest that the density waves (traveling density fiuctuations) are kinetic in nature. 
Furthermore, the fiux-density curve for the granular fiow resembles that of a traffic fiow, 
which is considered as a prime example of the systems which shows kinetic waves 0. 

One additional piece of evidence that the density waves are of a kinetic nature is their 
dependence on the initial density po- The theory of kinetic waves predicts 0] that small 
density fiuctuations in a uniform density background Po travel with a velocity 



which is the slope of the fiux-density curve at the mean density. We thus expect a 
large negative velocity for small po, a decrease to zero velocity at po ~ 15, with an in- 
creasingly large positive velocity as po is increased further. To check this, we measure 
the wave velocities for several values of po (keeping all other parameters fixed as above). 
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Writing the mean density po and the corresponding velocity U{po) as {po,U{po)), we find 
(10.0, -41 ± 2), (15.0, 5 ± 9), (18.7, 12 ± 11) and (22.5, 113.0 ± 4), which are all consistent 
with the above prediction. 

We now consider the theoretical aspect of the density waves. Consider the equations of 
motion which govern the time evolution of the density p{x,t) and the velocity v{x,t) fields 
for a granular fiow in a vertical tube. The first equation is that of mass conservation 

|p+^(P-)=0, (3a) 

and the second is a momentum conservation equation 

d d 

P^^ + P^^^ = ^(a;,t), (3b) 

where F{x,t)dx is the total amount of force acting on the particles in a region [x,x + dx]. 
The force F{x, t) has three contributions — gravity, internal pressure, and friction from the 
sidewalls. The exact form of the internal pressure and the friction is not known. Here, we 
use Bagnold's law, which is believed to be valid in the grain inertia regime Therefore, 
the force F{x, t) is 

d 

F{x, t) = -pgW - sign{v) pBfxyipy - D — [pBfxMv\ (4) 

Here, g is gravitational acceleration, ps the density of the material which forms the particles, 
p the packing fraction (p = Pbp) and D is the diameter of the particles. We assume the 
thickness of the shear layer to be of order of D. Also, f^x and f^y are geometry dependent 
functions, which contain the information about the density dependence of the forces. 
The uniform density solution of Eq. (|^) with the force given by Eq. (^) is 

p(x, t) = pbPo 

V{x,t) = - ^PogW/ f^y{Po) (5) 

If we add small density fiuctuation p = po + dp m the uniform density fiow, it can be shown 
|T0| that the fiuctuation travels with a velocity 



^ Jxy[Po) 

Equations and (|^) are exactly what one expects if the kinetic wave theory is to apply - 
uniform flow is a solution to the equations of motion, and density fluctuations travel with a 
density dependent velocity. 

Thus, it is clear then that the motion of the density pattern can be understood by 
applying the ideas of kinetic wave theory. However, this basic formalism only describes the 
motion of a pre-existing density pattern. It does not explain the observation that regions with 
large density contrasts are being formed out of the uniform background. Our simulations 
show that the large scale density pattern begins as a collection of small fluctuations in the 
density. These small fluctuations grow in time and a pattern emerges in which large density 
contrasts exist between neighboring regions. The evolution to such a state can be understood 
by considering the system as set of interacting kinetic waves. A detailed treatment of the 



general problem of interacting kinetic waves can be found elsewhere ||TT|]; in this paper, we 
present only the results from a simple model for the pattern evolution process in sand flowing 
down a tube. 

Consider the early stages of the flow in which the density of sand is nearly uniform at 
p ~ Po- Because of the roughness of the grains, the roughness of the walls or from the 
stochastic nature of the inelastic collisions, small density fluctuations appear in the system. 
In the interacting density wave approach, we treat the fluctuations as a set of distinct density 
regions with interfaces whose velocities are determined by a discrete form of Eq. (Q). In this 
case, the interface separating a region of density pi from a region with density p2 moves 
with a velocity, f/(l,2), given by 

t/(l,2) = ^-^^il^^, (7) 

Pi - P2 

which is the kinetic wave theory result for interfacial velocities involving finite density differ- 
ences The evolution of the system is determined by the motion of the interfaces and, as 
shall be shown, the nature of their interactions leads to a final state in which large density 
contrasts occur. 



In the computational and analytic results that follow, we choose a specific form for the 
density fiuctuations in the system. In our model, the initial positions of the interfaces are 
taken as a set of Nq points placed randomly on the interval [0,L], with regions between 
successive interfaces being assigned a density randomly in the range [po — W, po + W] . The 
principal virtues of this model are its simplicity and the fact that there are no correlations 
in the initial state which might infiuence the final structure. A more realistic model for the 
fiuctuations of the system would require a microscopic understanding of each specific source 
of noise. 

It is also necessary to choose a form for the fiux curve j (p). We have taken the parabolic 
form 



where R is the density at which no fiow occurs, and Jo is one quarter of the maximum fiux 
of the system. This curve was chosen for several reason. The first is that its simplicity 
eases some of the hardships of analytical calculations. The second reason is that for density 
fiuctuations over a sufficiently small range, the true fiux response can be approximated by 
this form (with R and Jo being fitting parameters). And finally, it is a first approximation 
to the form observed for the j{p) observed in Fig. 

Numerical simulation of this system is a very straight-forward exercise. The values of 
the densities in two adjacent regions determine the velocity of the associated interface. 
Consider three successive density regions A, B and C. During the course of the simulation, 
the interface A-B may encounter the interface B-C. This indicates that all of the mass that 
was inside region B has been completely "swallowed up" by the regions A and C. In this 
case, the interfaces A-B and B-C are replaced by a single A-C interface. The velocity of this 
interface can be calculated from the densities in regions A and C. Thus, it is a matter of 
tracking all of the interfaces, checking for collisions, and when they occur, replacing the two 
old interface with a single new one. Therefore, this technique does not allow for any density 
values other than those initially present, and the number of interfaces is always decreasing. 




(8) 
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For convenience, the simulations were done using periodic boundary condition. 

The first set of results shown below are from a simulation in which there are 400 interfaces 
initially placed randomly in the interval [0, 1] (i.e. L = 1). We also choose the values Jo = 1 
and R = 1. The densities are chosen at random from the interval [0.3, 0.8] (i.e. po = .55R, 
W = .25R). Figure ^(a) shows the initial density configuration, while figure |^(b) shows the 
system after a time t = .486 (where time is measured in the units of RL/Jo), and there are 
only 33 interfaces which remain along the interval. The system has evolved to a state in 
which the density contrast is very high between neighboring regions, and this behavior was 
observed for all values of po and W. 

This increase in the density contrast can be characterized quantitatively in the following 
way. Let the density of each region be pi, with i indexing the different regions, and N{t) be 
the number of regions at time t. Define the quantity 

^W = ]^ElP^-Pml, (9) 

where PN{t)+i = Pi- The larger the value of M(t) the larger the density contrast between 
neighboring regions. Figure ^ shows the quantity M{t) — M(0) averaged over 10 simulations 
with A'o = 10,000 interfaces. At early times, there is a linear increase in M{t) with a 
crossover to a nearly constant value at late times. 

At early times, it is possible to calculate M{t) analytically and the results are shown 
as the dotted line in Figure ^. In this regime, the changes in M{t) are dominated by the 
interaction of interfaces whose movements are determined by the initial configuration of 
the system. The calculation averages over all possible configurations of the initial random 
densities and interfaces, determines the time at which each interface collision occurs and 
how much that collision changes the value of M{t). In this regime, the agreement with 
the simulation is good. It is also possible to show exactly that M(0) = 2W/3. At later 
times, after there have been many collisions between interfaces, the structure of the system 
depends on the nature of the earlier evolution. Thus, this long time behavior is much more 
difficult to calculate. The results from the calculation described above break down in this 
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regime because the distribution of density regions is no longer that of the initial random 
distribution. 

At long times M{t) ^ 2W. Thus, the density contrast at long times is, on average, 
as large as the largest density contrasts present in the initial configuration. It turns out 
that the interacting kinetic waves do not create large contrasts. Rather, the interfaces from 



the initial distribution which survive are those that have a very large density contrast [|Tl 
Thus, while the noise in the system may provide a variety of such contrasts, the interacting 
kinetic waves will destroy all but the very largest. 

This paper outlines a kinetic wave approach to understanding the density patterns ob- 
served in sand flow along a vertical tube (many of the details omitted here can be found 



in references [|T0[ and ||TT|). However, these ideas certainly do not constitute a complete 
theory for the patterns observed in the experimental system. The role that the flow of air 
plays in this process |]12[, as well as the sources of noise in the system, are certainly not well 



understood. Further experimental investigation of these issues would be most enlightening. 
From a theoretical point of view, it is not clear whether the frictional force at the wall and 
the internal pressure obey Bagnold's law. While this form has been observed in the sheer 
cell geometry [0,0 , there has been no direct measurement of the frictional force for gravity 
driven flow. Finally, it is known that the interface between two regions of differing densities 
may not be a stable structure 0, and that diffusive effects may strongly influence the long 
time behavior of a system of interacting kinetic waves. 

The authors would like to thank the members of the HLRZ Many Body Group for 
stimulating discussion throughout this work. 
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FIGURES 

FIG. 1. Time evolution of (a) density and (b) velocity fields. These simulations were done with 
a tube of width W = 1 and length L = 15. Fields at a given time are shown as a vertical line 
of small boxes. The grayscale of each box is proportional to the value of the density or velocity 
in that region of the tube. Regions of high density are formed, and travel with almost constant 
velocity. 

FIG. 2. Local flux as a function of local particle density. This curve was found for a tube 
with width W = 1 and length L = 15, obtained by time averaging. The parabolic shaped curve 
resembles the flux-density relation in traffic flows. 

FIG. 3. Evolution of interacting kinetic waves. Both strips use a linear grayscale with white 
indicating po = and black po = R, the jamming density, (a) Shows the initial configuration of 
400 interfaces with densities in the range [0.3R,0.8R]. (b) Shows the configuration when only 33 
interfaces remain, and illustrates the tendency for alternating high and low density regions. 

FIG. 4. Contrast, M{t), as function of time, t. The solid line shows the results from averaging 
over 10 simulations with iV(0) = 10,000, and densities chosen in the range [0.3, 0.8]. At early times, 
the increase in contrast is linear and at long times it becomes a constant. The dashed line shows 
the results from an analytical calculation of the short time behavior. 
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